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A stochastic background of gravitational waves is expected to arise from a 
superposition of a large number of unresolved gravitational-wave sources of as- 
trophysical and cosmological origin. It is expected to carry unique signatures 
from the earliest epochs in the evolution of the universe, inaccessible to the 
standard astrophysical observations [1]. Direct measurements of the amplitude 
of this background therefore are of fundamental importance for understanding 
the evolution of the universe when it was younger than one minute. Here we 
report direct limits on the amplitude of the stochastic gravitational-wave back- 
ground using the data from a two-year science run of the Laser Interferometer 
Gravitational- wave Observatory (LIGO) |2Q. Our result constrains the energy 
density of the stochastic gravitational-wave background normalized by the crit- 
ical energy density of the universe, in the frequency band around 100 Hz, to be 
less than 6.9 x 10 6 at 95% confidence. The data rule out models of early universe 
evolution with relatively large equation-of-state parameter |3], as well as cosmic 
(super) string models with relatively small string tension [4j] that are favoured 
in some string theory models [5j. This search for the stochastic background 
improves upon the indirect limits from the Big Bang Nucleosynthesis jl[lg] and 
cosmic microwave background |7j] at 100 Hz. 

According to the general theory of relativity, gravitational waves (GWs) are produced 
by accelerating mass distributions with a quadrupole (or higher) moment. Moreover, in the 
early phases of the evolution of the universe, they can be produced by the mechanism of 
amplification of vacuum fluctuations. Once produced, GWs travel through space-time at the 
speed of light, and are essentially unaffected by the matter they encounter. As a result, GWs 
emitted shortly after the Big Bang (and observed today) would carry unaltered information 
about the physical processes that generated them. These waves are expected to be generated 
by a large number of unresolved sources, forming a stochastic gravitational- wave background 
(SGWB) that is usually described in terms of the GW spectrum: 

where dpcw is the energy density of gravitational radiation contained in the frequency 
range / to f + df and p c is the critical energy density of the universe [8|. Many cosmological 
mechanisms for generation of the SGWB exist, such as the inflationary models , pre- 
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FIG. 1: Sensitivities of LIGO interferometers. LIGO interferometers reached their design 
sensitivity in November 2005, resulting in the interferometer strain noise at the level of 3 x 10~ 22 
rms in a 100 Hz band around 100 Hz. This figure shows typical strain sensitivities of LIGO interfer- 
ometers during the subsequent science run S5. Also shown is the strain amplitude corresponding 
to the upper limit on the GW energy density presented in this paper (gray dashed line). Note 
that this upper limit is ~ 100 times lower than the individual interferometer sensitivities, which 
illustrates the advantage of using the cross-correlation technique in this analysis. 

big-bang models ll|, [l2|, ll3j , electroweak phase transition [141 ]. and cosmic strings 

16]. There are also astrophysical mechanisms, such as due to magnetars 1^] or rotating 



neutron stars 



The physical manifestation of GWs consists of stretching and compressing the spatial 
dimensions orthogonal to the direction of wave propagation, producing strain in an oscillating 
quadrupolar pattern. A Michelson interferometer with suspended mirrors Q is well suited 
to measure this differential strain signal due to GWs. Over the past decade, LIGO has built 
three such multi-kilometer interferometers, at two locations |2j: HI (4 km) and H2 (2 km) 
share the same facility at Hanford, WA, and LI (4 km) is located in Livingston Parish, LA. 
LIGO, together with the 3 km interferometer Virgo 19J in Italy and GEO 20j] in Germany, 
forms a network of GW observatories. LIGO has completed the science run S5 (between 
November 5, 2005 and September 30, 2007), acquiring one year of data coincident among 
HI, H2 and LI, at the interferometer design sensitivities (Fig. [1]). 

The search for the SGWB using LIGO data is performed by cross- correlating strain data 
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from pairs of interferometers [8j. In the frequency domain, the cross-correlation between 
two interferometers is multiplied by a filter function Q(f) (c.f. Data Analysis Supplement): 

This filter optimizes the signal-to-noise ratio, enhancing the frequencies at which the 
signal of the template spectrum ^gw(/) is strong, while suppressing the frequencies at 
which the detector noise {Pi{f) and i"2(/)) is large. In Eq. [2} and throughout this letter, we 
assume the present value of the Hubble parameter H = 72 km/s/Mpc [2l|, and use 7(/) 
to denote the overlap reduction function [8|, arising from the overlap of antenna patterns 
of interferometers at different locations and with different orientations. For the Hl-Ll and 
H2-L1 pairs the sensitivity above roughly 50 Hz is attenuated due to the overlap reduction. 
Since most theoretical models in the LIGO frequency band are characterized by a power 
law spectrum, we assume a power law template GW spectrum with index a: Qqw(/) = 
Q a (f/100 Hz) a . The normalization constant M in Eq. [2] is chosen such that the expected 
value of the optimally filtered cross-correlation is Q a . 

We apply the above search technique to the data acquired by LIGO during the science run 
S5. We include two interferometer pairs: Hl-Ll and H2-L1. Summing up the contributions 
to the cross-correlation in the frequency band 41.5-169.25 Hz, which contains 99% of the 
sensitivity, leads to the final point estimate for the frequency independent GW spectrum 
(a = 0): Q = (2.1 ± 2.7) x 10~ 6 , where the quoted error is statistical. We calculate the 
Bayesian 95% confidence upper limit for Q , using the previous LIGO result (S4 run [221 ]) as 
a prior for Q and averaging over the interferometer calibration uncertainty. This procedure 
yields the 95% confidence upper limit Qq < 6.9 x 10~ 6 . For other values of the power index a 
in the range between —3 and 3, the 95% upper limit varies between 1.9 x 10~ 6 and 7.1 x 10 -6 . 
These results constitute more than an order of magnitude improvement over the previous 
LIGO result in this frequency region 22j. Fig. [2] shows this result in comparison with other 
observational constraints and some of the cosmological SGWB models. 

Prior to the result described here, the most constraining bounds on the SGWB in the 
frequency band around 100 Hz came from the Big-Bang-Nucleosynthesis (BBN) and from 
cosmic microwave background (CMB) measurements. The BBN bound is derived from the 
fact that a large GW energy density at the time of BBN would alter the abundances of the 
light nuclei produced in the process. Hence, the BBN model and observations constrain the 
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FIG. 2: Comparison of different SGWB measurements and models. The 95% upper limit 
presented here, < 6.9 x 10 -6 (LIGO S5), applies in the frequency band 41.5-169.25 Hz, and is 



compared to the previous LIGO S4 result |22|] and to the projected Advanced LIGO sensitivity 
Note that the corresponding S5 95% upper bound on the total gravitational-wave energy density 
in this band, assuming frequency independent spectrum, is 9.7 x 10 -6 . The indirect bound due to 
BBN [l, 6] applies to ^bbn = / ^Gw(/)^( m /) (and not to the density ^gw(/)) over the frequency 
band denoted by the corresponding horizontal line, as defined in Equation A similar integral 
bound (over the range 1CP 15 - 10 10 Hz) can be placed using CMB and matter power spectra Q]. 
rojected sensitivities of the satellite-based Planck CMB experiment [7J and LISA GW detector 



261 ] are also shown. The pulsar bound 



27J is based on the fluctuations in the pulse arrival times of 
millisecond pulsars and applies at frequencies around 10~ 8 Hz. Measurements of the CMB at large 
angular scales constrain the possible redshift of CMB photons due to the SGWB, and therefore 



limit the amplitude of the SGWB at largest 



inflationary 
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. Examples of 



13 ] models are also shown 



161 ] , and pre-big-bang 

(the amplitude and the spectral shape in these models can vary significantly as a function of model 
parameters). 



total GW energy density at the time of nucleosynthesis 



n 



BBN 



J fi GW (/) d(\nf) < 1.1 x 1(T 5 (N v - 3), 



(3) 



where N u (the effective number of neutrino species at the time of BBN) captures the uncer- 
tainty in the radiation content during BBN. Measurements of the light-element abundances, 



11 




BBN (£2= 1.5x|10" = ) 
-LIGO (f2 = 6.9>|l0" 6 ) 
-AdvLIGO (n = \.0x 10~ 9 

Pulsar (i2 = 3.9|x 10 

i i i 

0.2 0.4 



0.6 



0.8 



FIG. 3: Constraining early universe evolution. The GW spectrum ^gw(/) is related to the 
parameters that govern the evolution of the universe 3J: 

Ogw(/) = ^4 /"W /"■*(/) r, where a(f) = 2 , r is the ratio of tensor and scalar perturbation 

amplitudes (measured by the cosmic microwave background (CMB) experiments), fi-t(f) and w(f) 
are effective (average) tensor tilt and equation of state parameters respectively, and A is a constant 
depending on various cosmological parameters. Hence, the measurements of ficw an d r can be 
used to place constraints in the w — fit plane, independently of the cosmological model. The figure 
shows the w — fit plane for r = 0.1. The regions excluded by the BBN 23J, LIGO, and pulsar 



bounds are above the corresponding curves (the inset shows a zoom-in on the central part of the 

3]. We note that the CMB bound [7] almost exactly 



figure). The BBN curve was calculated in 
overlaps with the BBN bound. Also shown is the expected reach of Advanced LIGO 



25|]. Note 



that these bounds apply to different frequency bands, so their direct comparison is meaningful only 
if fit(f) and w(f) are frequency independent. We note that for the simplest single-field inflationary 
model that still agrees with the cosmological data, with potential V{4>) = m 2 (j) 2 /2 (where <fi is a 



scalar field of mass m), r = 0.14 and nf(100 Hz) 
equation-of-state parameter of w(100 Hz) < 0.59. 



-0.035 



281 ]. implying a LIGO bound on the 



combined with the Wilkinson Microwave Anistropy Probe (WMAP) data give the upper 



bound N„ - 3 < 1.4 



231 ] . Similarly, a large GW background at the time of decoupling of 



CMB would alter the observed CMB and matter power spectra. Assuming homogeneous 
initial conditions, the total GW energy density at the time of CMB decoupling is constrained 
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to / Qgw(/) d(]nf) < 1.3 x 10" 5 Q. In the LIGO frequency band and for a = 0, these 
bounds become: Qq BN < 1.1 x 10~ 5 and f2Q MB < 9.5 x 10~ 6 . Our result has now surpassed 
these bounds, which is one of the major milestones that LIGO was designed to achieve. 
Moreover, the BBN and CMB bounds apply only to backgrounds generated prior to the 
BBN and the CMB decoupling respectively, while the LIGO bound also probes the SGWB 
produced later (this is the case, for example, in cosmic strings models). 

Our result also constrains models of the early universe evolution. While the evolution 
of the universe following the BBN is well understood, there is little observational data 
probing the evolution prior to BBN, when the universe was less than one minute old. The 
GW spectrum flcwif) carries information about exactly this epoch in the evolution. In 
particular, measuring G w(/) is the best way to test for existence of presently unknown 
"stiff" energy components in the early universe jjj, for which a small density variation is 
associated with a large pressure change, which could carry information about the physics of 
the inflationary era [24j . Fig. [3] demonstrates how the result presented here can be used to 
constrain the existence of these new energy components. 

Our result also constrains models of cosmic (super) strings. Cosmic strings were origi nally 
proposed as topological defects formed during phase transitions in the early universe [15] . 
More recently, it was realized that fundamental strings may also be expanded to cosmological 
scales j^j. Hence, searching for cosmic strings may provide a unique and powerful window 
into string theory and into particle physics at the highest energy scales. Fig. H] shows 
that our result, along with other observations, can be used to constrain the parameters in 
the cosmic string models. While our result is currently excluding a fraction of the allowed 
parameter space, Advanced LIGO [25| is expected to probe most of these models. 

Measurements of the SGWB also offer the possibility of probing alternative models of the 



early universe cosmology. For example, in the pre-Big-Bang model [ll|, [121, [13j the universe 
starts off large and then undergoes a period of inflation driven by the kinetic energy of a 
dilaton field, after which the standard cosmology follows. Although more speculative than 
the standard cosmology model, the pre-Big-Bang model makes testable predictions of the 
GW spectrum. As shown in Fig. [51 the BBN and CMB bounds are currently the most 
constraining for this model and Advanced LIGO {25] is expected to surpass them. 
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FIG. 4: Cosmic strings models. The network of cosmic strings is usually parametrized by the 
string tension fx (multiplied by the Newton constant G), and reconnection probability p. The 
CMB observations limit Gji < 10~ 6 . If the size of the cosmic string loops is determined by the 
gravitational back-reaction 



161 ] which 



29l | , the size of the loop can be parametrized by a parameter e 
is essentially unconstrained. The mechanism for production of GWs relies on cosmic string cusps: 
regions of string that move at speeds close to the speed of light. If the cusp motion points toward 

|3Q |. The superposition 



Earl 1 1. a dried able burst, of gravitational radiation may be produced [lfl . 
of GWs from all string cusps in the cosmic string network would produce a SGWB 4|. This 
igure shows how different experiments probe the e — G/i plane for a typical value of p = 10 -3 
4| (p is expected to be in the range 10 -4 — 1). The excluded regions (always to the right of the 



corresponding curves) correspond to the S4 LIGO result [22]. this result. IJBN bound 



6, 



2J, CMB 



bound [7|, and the pulsar limit [27]]. In particular, the bound presented in this paper excludes a 
new region in this plane (7 x 10~ 9 < G/i < 1.5 x 1CP 7 and e < 8 x 10 ), which is not accessible to 
any of the other measurements. Also shown is the expected sensitivity for the search for individual 
bursts from cosmic string cusps with LIGO S5 data [301 ] . The region to the right of this curve 
is expected to produce at least one cosmic string burst event detectable by LIGO during the S5 
run. Note that this search is complementary to the search for the SGWB as it probes a different 
part of the parameter space. Also shown is the region that will be probed by the Planck satellite 
measurements of the CMB {tJ. The entire plane shown here will be accessible to Advanced LIGO 



25[ SGWB search. 
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FIG. 5: Pre-Big-Bang models. In the pre-Big-Bang model, the GWs are produced via the 
mechanism of amplification of vacuum fluctuations, analogously to the standard inflationary model. 
The typical GW spectrum increases as / 3 up to a turn-over frequency f s , above which Hgw(/) ~ 
y3-2^t w ^h n < 1.5. The spectrum cuts off at a frequency fi, which is theoretically expected to be 
within a factor of 10 from 4.3 x 10 10 Hz (dashed horizontal line). This figure shows the f\ — fi plane 
for a representative value of f s = 30 Hz. Excluded regions corresponding to the S4 result and to 
the ^ rented here » Shaded. The r e gi „»s exceed h y the BBN Q,Q and the CMB Q 
bounds are above the corresponding curves. The expected reaches of the Advanced LIGO [25] and 
of the Planck satellite [7] are also shown. 



I. DATA ANALYSIS SUPPLEMENT 



A. Method 



Gravitational waves stretch and compress the spatial dimensions perpendicular to the 
direction of wave propagation. In a Michelson interferometer with suspended mirrors, the 
gravitational wave would cause stretching and shrinking of orthogonal arms, as shown in 
Figure El which would result in corresponding fluctuations in the laser intensity at the 
output of the interferometer. Hence, transient or periodic gravitational waves would cause 
transient or periodic fluctuations in the output laser power. A stochastic gravitational- 
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FIG. 6: Effect of a gravitational wave on an interferometer. A gravitational wave traveling 
perpendicular to the plane of the interferometer stretches and compresses interferometer arms in an 
alternating manner. The laser beam (entering from the left) is split equally between the arms, the 
two new beams travel to and reflect back from the end mirrors, and are superposed at the photo- 
detector (at the bottom). Changes in the arm lengths cause the two beams to acquire different 
phases while traveling in the arms, the differential component of which is observed as modulations 
in the laser light intensity at the photo-detector. 

wave background (SGWB) signal would cause random fluctuations in output laser power, 
which are indistiguishable from various instrumental noise sources. We hence search for a 
SGWB by cross-correlating strain data from pairs of interferometers, as described in [8]. In 
particular, we define the following cross-correlation estimator: 



where 5t is a finite-time approximation to the Dirac delta function, si and §2 are the 
Fourier transforms of the strain time-series of two interferometers, and Q is a filter function. 
Assuming that the detector noise is Gaussian, stationary, uncorrelated between the two 
interferometers, and much larger than the GW signal, the variance of the estimator Y is 
given by: 



where Pi(f) are the one-sided strain power spectral densities (PSDs) of the two interferom- 
eters and T is the measurement time. Optimization of the signal-to-noise ratio leads to the 




(4) 
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following form of the optimal filter [8]: 



W)=AT 



7(/)fW/)# 2 
/ 3 Pi(/)P 2 (/) 



(6) 



where H is the present value of the Hubble parameter, assumed below to take the value H = 
72 km/s/Mpc [19], and 7(/) is the overlap reduction function [8], arising from the overlap of 
antenna patterns of interferometers at different locations and with different orientations. For 
the Hanford-Livingston pairs the sensitivity above 50 Hz is attenuated due to the overlap 
reduction, while the identical antenna patterns of the colocated Hanford interferometers 
imply 7(/) = 1. Hence, the colocated Hanford interferometer pair is more sensitive to 
the isotropic SGWB than the Hanford-Livingston pairs, but it is also more susceptible to 
environmental and instrumental correlations. For this reason, this pair is not included in 
the analysis presented here. Since most theoretical models in the LIGO frequency band are 
characterized by a power law spectrum, we assume a power law template GW spectrum 
with index a, 



The normalization constant Af in Equation [6] is chosen such that < Y >= Q a . 
B. Results 

Our results are based on the LIGO data acquired during the science run S5, which 
took place between November 5, 2005 and September 30, 2007. Virgo [19] and GEO [20] 
detectors were also operating during some parts of this science run. However, due to their 
lower strain sensitivities around 100 Hz, these interferometers were not included in the 
search presented here. We analyzed the Hl-Ll and H2-L1 interferometer pairs. The data 
for each interferometer pair was divided into 60 sec segments, down-sampled to 1024 Hz, 
and high-pass filtered with a 6 th order Butterworth filter with 32 Hz knee frequency. Each 
segment I was Hann-windowed and estimators Yj(f) and cri(f) were evaluated with 0.25 Hz 
resolution. To recover the loss of signal-to-noise due to Hann- windowing, segments were 50% 
overlapped. A weighed average was performed over all segments from both interferometer 
pairs, with inverse variances as weights. 

The data were preselected to avoid digitizer saturation effects, periods with unreliable 
calibration, and periods suffering from known instrumental transient disturbances. In ad- 
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FIG. 7: Histogram of the fluctuations of the estimator Yj over segments / around the mean, 
normalized by the standard deviation 07 is shown in red (for the Hl-Ll pair). The blue curve 
shows the Gaussian fit to the histogram, which has zero mean and unit variance. The Kolmogorov- 
Smirnov test statistic (comparing the histogram and the fit) is 0.2 for Hl-Ll (0.4 for H2-L1), 
indicating that the data is indeed Gaussian-distributed, and that the estimate of the theoretical 
variance a\ is reliable. 

dition, about 3% of the segments were found to deviate from the assumption of stationary 
noise: the difference between 07 and a calculated using the neighboring segments exceeded 
20% for these segments, and they were not included in the analysis. The 20% threshold is 
optimal as it yields gaussian distribution of the data (c.f. Figure [7]), while minimizing the 
amount of eliminated data. The data quality selection was performed blindly, using an un- 
physical 0.5-sec time-shift between the two interferometers (a broadband SGWB covering 
the range of ~ 100 Hz is expected to have coherence time ~ 10 ms, as also depicted in Figure 
ITU1) . Once the data selection was completed, the final zero-lag analysis was performed. The 
selected segments amount to 292 days of exposure time for Hl-Ll (294 days for H2-L1). 

To identify potentially contaminated frequency bins, we calculated the coherence between 
HI and LI (and H2 and LI) over the entire S5 run. The coherence is defined as 



where (Pn(f)) is the average strain cross-spectral density between two interferometers and 



r(/) 




(8) 
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FIG. 8: Coherence between HI and LI strain data is shown at 1 mHz resolution (top) and 100 mHz 
resolution (bottom). The insets show the histograms of the coherence along with the expected 
exponential distribution. Note that after notching the contaminated bins (red), the remaining 
frequencies follow the expected exponential distribution. Note: N denotes the number of averages 
used in the calculation. 

(Pi(f)) is the average strain power-spectral density for the interferometer %. These calcula- 
tions have revealed several instrumentally correlated lines between each pair of interferom- 
eters: 16 Hz harmonics (associated with the data acquisition clock), 60 Hz harmonics (AC 
power line), and injected simulated pulsar signals (52.75 Hz, 108.75 Hz, 148 Hz, 193.5 Hz, 
and 265.5 Hz). These lines were found to be correlated between instruments in the blind 
analysis, and were excluded from the final zero-lag analysis. Figure [8] shows the coherence 
between HI and LI strain data at 1 mHz and 100 mHz resolutions. 
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FIG. 9: Stochastic signal simulations in hardware for Hl-Ll (blue) and H2-L1 (red), and in software 
(Hl-Ll, green) are shown. The error bars denote 2a ranges. 
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FIG. 10: Signal-to-noise ratio for the recovery of a software simulation with Hl-Ll data with 
(^simulated _ 3 g x j g s j lown as a function of the time-lag between the two interferometers. 
The inset shows the zoom-in around zero-lag: the signal is recovered well for zero-lag (SNR ~ 7.2), 
but it disappears quickly with time-lag of ±30 ms. 

The search algorithm described here is verified using signal simulations. The simulations 
are performed in hardware (by physically moving the interferometer mirrors coherently be- 
tween interferometers), in which case they are short in duration and strong in amplitude. 
They are also performed in software, by adding the stochastic signal to the interferometer 
data, in which case they can be long in duration and relatively weak in amplitude. Three 
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FIG. 11: Y(f) and tr(/) obtained by combining the Hl-Ll and H2-L1 data from the S5 run. The 
inset shows that the ratio of the two spectra is consistent with a Gaussian of zero mean and unit 



variance. 
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FIG. 12: Upper limit is shown as a function of the power index a for several LIGO results: based 
on the previous runs S3 and S4 and the S5 result presented here. 

hardware simulations were performed, with amplitudes of Qq ~ 2 (20 min long), 2 x 10 -2 
(20 min long), and 6.5 x 10~ 3 (~ 3.8 hours long) and they were successfully recovered (within 
experimental uncertainties) for both Hl-Ll and H2-L1 pairs. A software simulation was per- 
formed and successfully recovered using about 1/2 of the Hl-Ll data, with the amplitude of 
Qo = 3.8 x 10~ 5 . Figures 191 and [TUl demonstrate the recovery of both hardware and software 
simulations. 
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We apply the above search technique to the data acquired by LIGO during the science 
run S5, repeating it for the interferometer pairs Hl-Ll and H2-L1. We treat the data from 
the two pairs as uncorrelated, although HI and H2 are known to suffer from instrumental 
and environmental correlations. We have verified that the level of the H1-H2 correlations is 
sufficiently small that it could affect the result presented here by less than 1%. The resulting 
composite spectrum for the frequency independent template (a = 0) is shown in Figure [Til 
Integrated over the frequency band 41.5-169.25 Hz, which contains 99% of the sensitivity as 
determined by the variance integrand, this leads to the final point estimate for the frequency 
independent GW spectrum: Q = (2.1 ± 2.7) x 1CT 6 , where the quoted error is statistical. 
We calculate the Bayesian posterior distribution for Q using this result. For the prior 
distribution of Qq we use our previously published posterior distribution from the earlier 
S4 run [22]. We also marginalize over the calibration uncertainty, which is the dominant 
systematic error in this analysis and was estimated to be 13.4% for LI and 10.3% for HI 
and H2. With these assumptions, the final 95% confidence upper limit is Qo < 6.9 x 1CT 6 . 
Figure [12] shows the 95% confidence upper limit as a function of the power index a of the 
template spectrum. This result is more than an order of magnitude improvement over the 
latest LIGO result in this frequency region [22]. 

C. Outlook 

LIGO and Virgo are planning significant upgrades to their interferometers, known as 
Advanced LIGO and Advanced Virgo. These upgrades will improve the interferometers' 
strain sensitivities by 10 times across the entire frequency band, and they will extend the 
sensitive frequency band down to ~ 10 Hz. Consequently, the network of advanced detectors 
will be able to probe the isotropic SGWB at the level of flaw ~ 10~ 9 or better. Moreover, 
while searches for isotropic SGWB tend to be dominated by pairs of nearby or co-located 
detectors, the presence of the third location in the network is crucial for searches for non- 
isotropic SGWB. Hence, the network of advanced detectors is expected to produce detailed 
maps of the gravitational-wave sky, potentially revealing non-isotropic sources of SGWB, 
such as point sources or sources distributed in the galactic plane. Techniques for performing 
the searches for non-isotropic SGWB are currently under development. 
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